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Recent contributions in the field of quantum state tomography have shown that, despite the expo- 
nential growth of Hilbert space with the number of subsystems, tomography may still be performed 
efficiently by tailored reconstruction schemes. Here, we discuss a scalable method to reconstruct 
mixed states that are well approximated by matrix product operators. The reconstruction scheme 
only requires local information about the state, giving rise to a reconstruction technique that is 
scalable in the system size. It is based on a constructive proof that generic matrix operators are 
fully determined by their local reductions. 



The complexity of many-body systems is one of the 
most intriguing, but at the same time daunting, features 
of quantum mechanics. The "curse of dimensionality", 
namely the exponential growth of the descriptive com- 
plexity of even pure states, is a property of quantum 
mechanics which clearly distinguishes it from classical 
physics. Therefore, in general, the number of variables 
required to uniquely determine a quantum state increases 
in accordance with the growth of the Hilbert space expo- 
nentially. 

Quantum state tomography addresses the problem of 
completely characterizing a state of a physical system by 
measuring a complete set of observables that determine 
the state uniquely [1 . As the complexity of quantum op- 
erations implemented in a laboratory steadily increases 
[2]-[4] , the demand for a reliable and scalable tomography 
[5l[6] of prepared states is high and of considerable impor- 
tance for the future of quantum technologies. The ability 
to store and manipulate interacting quantum many-body 
systems, such as linearly arranged ions in an ion trap, 
enhanced rapidly during the last years such that large 
number of quantum systems will soon be controllable for 
a number of particles where the conventional methods of 
full quantum state tomography fail due to both, the lim- 
ited time that is realistically available for the experiment 
and the limitations to the resources that are available for 
the classical post-processing of the experimental data |2] . 

The experimental time requirement is defined by the 
amount of measurements which have to be done to recon- 
struct the state faithfully, i.e. one has to consider the sys- 
tem size and the number of repetitions to obtain sufficient 
statistics [7]. The latter is restricted by the individual 
tomography scheme and in particular by the representa- 
tion of the state. Clearly, full quantum state tomography 
where the state is represented by an exponentially large 
number of variables will require an exponentially increas- 
ing computational power and is hence infeasible already 
for, e.g., trapped ion experiments available today [4 . But 
naturally occurring quantum states and many states of 
interest for quantum information tasks are completely 



characterized by a number of variables scaling moder- 
ately in the number of particles: Ground states of gapped 
local Hamiltonians [SHTO] , thermal states of local Hamil- 
tonians [10^ JJJ, the W state, the GHZ state, and cluster 
states are all matrix product operators (matrix product 
states if they are pure) of low dimension, or very well ap- 
proximated by them. These states are parameterised by 
a linear number of matrices of low bond-dimension. The 
key insight here is not that these states are matrix prod- 
uct operators or states {any state is a matrix product 
operator, respectively state) but that the matrix dimen- 
sion is low, in particular independent of the system size. 
This solves one issue concerning the post-processing side 
of the problem mentioned above as these states may be 
stored efficiently on a classical computer. On the other 
hand, as we will see, generic matrix product operators are 
not only completely determined by a linear number of lo- 
cal observables but may also be efficiently reconstructed 
from such local measurements, which makes the formal- 
ism we present here a powerful tool for quantum state 
tomography. 

Recently, it has been demonstrated that the recon- 
struction of pure quantum states for large systems is 
possible with the knowledge of local information only [5] . 
Here, we develop new theoretical results that permit the 
generalization of these ideas to mixed states. We con- 
sider N (i-level subsystems aligned in a one-dimensional 
geometry, e.g., a chain of qubits {d = 2). The aim is 
to reconstruct a mixed state g from local information 
only. The local information we have in mind are esti- 
mates to all reductions of the state to a fixed number R 
of contiguous sites. These may be obtained by estimat- 
ing the expectation values of an informationally complete 
set of observables on the R sites. Note that we do not 
require the estimates to the reductions to be states, i.e., 
empirical estimates to the expectation values of local ob- 
servables suffice and post-processing such as maximum 
likelihood estimation is not required. As R is fixed, this 
corresponds to an experimental effort that is linear in 
the system size TV. We present a computationally cheap 
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tomography scheme, that is scahng polynomiahy in the 
system size which succeeds provably under a certain 
technical invertibility condition on g and demonstrate nu- 
mericahy that it stih works weh when this condition is 
not met. 



THE SCHEME 

We begin our exposition of the tomography scheme 
by introducing some notation. We denote the to-be- 
reconstructed state by g. The input to the reconstruction 
scheme are estimates of expectation values which com- 
pletely specify all reductions of ^ to contiguous sites. 
We denote these reductions by gk^ k = 1^ . . . ^ N — R-\- 1. 
Put mathematically, gk = tr{i,...,/,_i|u{fe+i?,...,iv} [^], i.e., 
the trace over all but the R sites {/c, . . . , + — 1}. 
Now let {P/"^}Q;=i,...,d2 be a complete operator basis for 
the site i. A common choice for spin- 1/2 particles is 



given by ^^'^ = 1,/V2, P^^'^ = af/V2, Pl'^ = /V2, 
/V^, i-e., the orthonormal Pauli spin basis. We 



may then write 



Qk 



\ p(«i) p("fi) (^\ 
k "'-^k+R-l-> 



i.e., the gk are completely specified by the local expecta- 
tion values (P^-) • • • Pili_,), = tr[^p(-) • • • P^^l J, 
estimates to which are the input to our tomography 
scheme. 

The gk completely specify the state ^ if a certain tech- 
nical invertibility condition is met. The proof is construc- 
tive and gives an explicit method for obtaining g from 
the gk. It is partly inspired by the characterization of 
finitely correlated states (as opposed to C*-finitely cor- 
related states) on infinite spin chains as described in the 
early literature [12, Prop. 2.1]. In addition to the fact 
that we are working in a finite and non-translation in- 
variant setting, the main novel technical point here is 
that we only use local information, provided by the gk. 
To state the invertibility condition, we first need to estab- 
lish some notation. We collect the TV sites of the system 
in the set A/" = {1, ... , N}. For X C A/", we define the 
complex vector spaces Vx spanned by 



= l,...,d2 



(2) 



For given O G V_\f and I, J C A/" we define the linear 
map -.Vj ^ V^\i as 



(3) 



We note that the map depends only on the reduction 
Oxuj = ^^Af\xuj[0] of O to sites XU J7 as 



this is illustrated in Fig. [T] Note that from now on we 
will only consider cases where XU J is connected. 

1 ••• k-l k k + r ••• N 

oo • ooooo oo, 

l-XH I J 1 

I Af 



FIG. 1. Definition of the sets A/" = {1, . . . , A^} and X, J' C A/". 
The linear map (X) = tij[XOiuj] maps operators X 
(e.g. observables) living on set J into operators on set X by 
means of the reductions of O (e.g. the state) to XU J'. 



Definition 1 (Invertibility) XeH, r G N, 2 < / + r < 
N-2. If 6 is such that for allk e¥lj<k< N-r-1, 
the equality 



rank 



E 



{/c+l,...,/c+r} 
{k-l-\-l,...,k} 



rank 



'{k+l,...,N} 



(5) 



holds, we call O (l^r)-invertihle. 



We may now state the main theorem, a proof of which 
may be found in the appendix. 

Theorem 1 Letl.r such that 2 <l^r < N - 2. 
Let O G Vj\f he {l,r)-invertible. Then, for all Xi G V^iy, 
the equality 



tTL[Xi---XNd]=tTL[Xi---XiYid] 



(6) 



holds. Here, the Yi G are recursively defined 

as follows. We set Yn-v = ^N-r-\-i ' • • ^^g? 



for k = / + 1, . . . , — r. Here, 
Moore-Penrose pseudoinverse. 



the bar indicates the 



Note that, for Eq. ([6|) the reduction of O to sites 
{1, + r} is needed, for the inverse we require the 
reduction of O to sites {/c — A: + r — 1}, and for 



{/c,...,/c+r} 
{k-l,...,k-l} 



^y{Xk^k) we require the reduction of O to sites 



E 

{/c — /,...,/c + r}. Hence, expectation values of the form 
tr/,[Xi • • • XatO] are completely specified by reductions 
of O to the sites {/c — /c + r}, /c = / + 1, . . . , TV — r, 

i.e., by all reductions toR = r-\-l-\-l contiguous sites. 
By choosing the Xi to be the basis operators P!^\ this 
implies that (/, r)-invertible operators O may be fully re- 
constructed from their reductions to R consecutive sites, 
which is the same as knowing the expectation values 



tr[4"^ 



k-\-R- 



(8) 



(4) for all /c = l,...,Ar-i?+ 1. 
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Now, of course, experimentally, the exact expectation 
values are only known to within a certain statistical er- 
ror (J (commonly the estimated standard deviation of 
the mean after a finite number of measurements). This 
error propagates into the singular values of the map 
E^^k-i^^k-l] ' ^^^^ needs to be inverted, even 
small errors on singular values close to zero will lead to 
a large error in the reconstruction. This issue may be 
avoided by using stochastic robust approximation tech- 
niques [15], see the appendix for details. 

In the appendix, we also show that a vast majority 
of matrix product operators fulfil the invertibility con- 
dition, i.e., a vast majority of matrix product operators 
may be reconstructed from local reductions alone. As 
noted above, practically relevant states are (well approx- 
imated by) matrix product operators of low dimension, 
i.e., we expect the scheme to work for a large class of 
mixed states. In the next section, we present numeri- 
cal results for states that do not necessarily fulfil the in- 
vertibility condition and for which the local expectation 
values are subject to inevitable statistical noise. 




FIG. 2. Quality of our reconstruction scheme for thermal 
states of the Ising Hamiltonian in Eq. ^ for /3 = 5 and 
i? = 5, i.e., the state is reconstructed from local expectation 
values on five consecutive sites. For each pair {N,a), the 
plot shows the mean of the norm difference obtained from 
100 realizations and renormalized by the purity of the exact 
state, i.e. D {g, Qrec) = \\grec - qW'^/WqW'^' This corresponds 
to 100 experiments, each of which carries an uncertainty of a 
about the local expectation values. 



NUMERICAL EXPERIMENTS 

We restrict our attention to qubits, d = 2, and illus- 
trate the behaviour of the tomography scheme for two 
different classes of mixed states: (i) Thermal states of 
local Hamiltonians and (ii) mixed states obtained by 
tracing out parts of a matrix product state. We sim- 
ulate the measurements in the following way. We first 
compute the exact local expectation values Pai,...,^^ ~ 



o^i = ^i^iUiZ^ for all k. Statisti- 



cal noise is then simulated by adding random numbers 
(drawn from a Gaussian distribution with zero mean and 
standard deviation a) to them. The resulting P^^^...^q,^ 
then serve as the input to our reconstruction scheme. 
At this point we would like to stress again that exper- 
imental estimates (here simulated by the ...an) 
the local expectation values suffice. Of course, these 
do not necessarily correspond to a state. We compare 
the reconstructed state ^^ec to the exact state g by com- 
puting the Hilbert-Schmidt norm difference D (^, ^rec) = 
ll^rec To obtain meaningful results, we have 

rescaled the norm such that the deviations are measured 
in units of the natural length scale of the state to 

be learned. 

As a first example, we consider thermal states of the 
Ising Hamiltonian at its quantum critical point 



H 



N-l 
i=l 



N 



i+1 



(9) 



We obtain the thermal states by an imaginary time 
evolution [13, 14 using the time evolving block decima- 
tion algorithm (TEBD). From these states, we compute 



the exact local expectation values p^^ 



for all /c, simulate the measurements by adding random 
numbers, and reconstruct the state by means of the noisy 
local expectation values. In Fig. [2j we show the norm 
difference for the exact and the reconstructed states as a 
function of the system size N and the error a. It indicates 
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FIG. 3. Quality of our reconstruction scheme for thermal 
states of randomly chosen next-neighbour Hamiltonians as in 
Eq. (10) with f3 — 2 and = 5, i.e. the state is reconstructed 
from local expectation values on five consecutive sites. Red: 
system size N = 16, blue: system size N = 32. We generate 
50 different random Hamiltonians and compute their corre- 
sponding thermal states using the TEBD algorithm. For each 
state and pair (A/", a), the density plot shows the simulations 
of one experiment carrying an uncertainty of a about the local 
expectation values. Mean values are indicated as triangles. 
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FIG. 4. Reconstruction errors for randomly chosen matrix product operators as described in the text with lA/luxl = N. The 
interaction is weak in a sense that we choose t such that t||^fc||op = 1/100 for all k. For each pair {N, a) we draw 4000 random 
states, simulate one measurement each and reconstruct the state with the disturbed local expectation values. The plot shows 
the mean values of the renormalized norm differences in dependence on the system size N and the error in the measurements a. 
(a) The states are reconstructed with = 3, i.e. measurements are done on all blocks of three contiguous sites. Here, for given 
N (cr), the scaling of D (g, Qrec) is roughly linear in a (N). (b) Reconstruction with = 5. D (g, Qrec) improves significantly 
when measuring on larger blocks. 



that, for given the error D ^rec) scales roughly as 
cr, similarly, for given a, it scales roughly as N. 

Thermal states of random next- neighbour Hamiltoni- 
ans of the form 

N-l 

serve as our second example. Here, the f • are Her- 
mitian matrices acting on sites i and i + 1 with entries 
that have real and imaginary part picked from a Gaus- 
sian distribution with zero mean and standard deviation 
one. Again, we use the TEBD algorithm to obtain the 
exact thermal states. For each system size we generate 
50 random Hamiltonians and their corresponding ther- 
mal states and simulate one experiment for each cr and 
state. Fig. |3] shows the error of the reconstructions as a 
function of the error of the measurements for two differ- 
ent system sizes. The densities illustrate the distribution 
of the error for the 50 different states while the black 
arrows indicate the mean. 

Let us finally analyse the behaviour of the algorithm 
for states which are exactly representable as matrix prod- 
uct operators satisfying the invertibility condition but 
subject to statistical noise. We pick such matrix prod- 
uct operators at random by generating a matrix product 
state with bond-dimension D = d where the entries of the 
matrices defining the states are drawn from a Gaussian 
distribution with zero mean and standard deviation one. 
Then, we let these sites interact with an auxiliary system 
each of dimension d according to the unitary Uk = e~*^^^ 
for A: = 1, ... , lATauxl, where Hk is a two-particle interac- 
tion Hamiltonian acting on site k and its auxiliary sys- 
tem with entries picked from a Gaussian distribution with 
zero mean and standard deviation one. Finally, we trace 



over the lA^auxl auxiliary sites to obtain a matrix prod- 
uct operator with bond-dimension D = d^ . Fig. [i] shows 
the results for different system sizes and different noise 
levels and further outlines the influence of enlarging the 
number of measured observables. Again, the numerical 
results suggest that the scaling of our scheme is polyno- 
mial in both, N and cr. 

CONCLUSIONS 

In this work we have continued the programme for effi- 
cient quantum state tomography initiated in [5] and pre- 
sented a scheme to reconstruct mixed states from local 
measurements efficiently. We have shown that practically 
relevant states may be reconstructed from reductions to 
contiguous sets of sites alone. Further, the impact of 
statistical noise and the performance of the reconstruc- 
tion scheme for states that do not necessarily fulfil the 
condition which guarantees perfect reconstruction have 
been investigated numerically in detail. For all simu- 
lations the Hilbert-Schmidt norm difference (normalized 
by the purity of the exact state) between the exact state 
and the reconstructed state was obtained and the numer- 
ical results suggest that the quality of the reconstruction 
scales algebraically in N and a. The methods presented 
here overcome an important challenge of tomography and 
open the path for their application to very large systems 
in concrete experiments. 
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APPENDIX 
Reconstructing invertible states 

Here, we proof the raain theorera. We start by showing 
that for ah = / + 1, . . . , A/" - r, Eq. ^ implies 



^{fe,...,/c+r-l}^ 



(rfc_i) = £;J^;;;,%>(x,n). 



(11) 



To this end, we define the linear map (j) : 
ran[£^l^'""''^^]^|] ran[^|^^y'^^^^2i].]? where the domain 
and the range of (j) are the ranges of the denoted linear 



maps, as 

p{fe,...,/c+r} 



7^{/e,...,/c+r} /r7\ 



(12) 



i.e., ran[^] = ran[£^|^^y'^^^^2i|]7 and therefore, by the 
rank-nullity theorem, 

dim[ker[</.]] = rank[^;,';-; - rank[sj,':y;'+;2,j] 



(13) 



due to the invertibility condition Eg . ([5| . Hence, = 
is equivalent to Z = 0, i.e., Eq. is equivalent to 



4>{e\Z:T.[\{Yu-i®^^^ (14) 

which is implied by Eq. The theorem now follows by 
induction over /c = A^ — r — 1,...,/ + 1, starting at /c = 
N — r — 1: Eq. ([s]) guarantees the existence of Yn-t-i ^ 
y{N-r,...,N-i} such that 

{N-r,...,N-l}(^ . _ j^{N-r,...,N} ^ x 



(15) 
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i.e., multiplying from the left by Xi • • • X^-r-i and tak- 
ing the trace over {1, . . . , — r — 1}, we find 

tr[Xi • • • XnO] = tr[Xi • • • X^nO] (16) 

for /c = A/" — r — 1. Suppose now that this equality holds 
for / < /c < A/" — r — 1 for some Yk ^ ^{/c+i,...,/c+r}- We 
now show that it then also holds for / < k — 1 < N — r — 2. 
Eq. ([s) guarantees the existence of Yk-i ^ ^{fe,...,/c+r-i} 
such that 



multiplying from the left by Xi • • • X/^-i and taking the 
trace over {1, . . . , /c — 1}, we find 



{/c,...,/c+r} / 



tr[Xi---X^O] =tr[Xi---MO] 

= tr[Xi---X,_in_iO], 

the desired equality for/</c — 1<A' — r — 2. 



(18) 



Generic matrix product operators are invertible 

Now consider matrix product operators 
0= A[c.i]---Piv[aiv]A^"^^---Plr"\ (19) 



with Pi[a] e Ci><^S PN[a] e C^^><\ and Pi[a] e 
(^DixDi+^ for z = 2, ... ,7V - 1. We assume w.l.o.g. that 
A^°^ (X 1, for ah z = l,...,Ar. 



Lemma 1 Let l,r e ¥1 such that 2</ + r<7V-2. 
Let O be a matrix product operator as in Eq. (19). If 
tr[d] and for all k e ¥1 , I < k < N - r - 1, the sets 



span (C^k-i+ixDk+i ^^g^ ^/^g 

{P/e+i[a/c+l] • • •P/c+r[Q^/c+r]}afc+i,...,afc+, (21) 

span (C^fc+ix^fc+r+i Qy^rp (j^^ then O is {I ^r) -invertible. 
Proof For X e V{/c+i,...,/c+r}, 

^ak+i,...,ak+r--^k+l '"-^k+r ' V^^/ 



we find 

^{/c+l,...,fe+r} 



'{/e-^ + l,...,/c} V"^/ 

^ A [1] • • • Pk-l [l]Pk-l + l Wk-l+l] ■■■Pk Wk] 



^afc + i,...,afc+r^fc+l[<^/c+l] • • ■ Pk+r[o^k+r] 



xp,+.+i[i].-.PM[i]prvjr---4^^ 

^ ii^'^'Pfe-^+i [ak-i^i] ■■■Pk [(^k]Xv 



Oik-l + l,---,Oik 



=:r(X), 



(23) 



where the matrix 

X = Xak+u...,ak+rPk+l[o^k+l] - - - Pk+r[o^k+r] 



ak+i,...,ak+r 



the vectors 

V = Pk+r+l[l]---PM[l] e C^-^+^+^Xl, 

«;t =Pi[l]---Pfc_,[l] eCi><^^-'+S 



(24) 



(25) 



and the mapping T : (D^fc+i x^^+^+i -> V{k-i+i,...,k}- 
Now, r(X) = is equivalent to 



= w'' Pk-i+i[ak-i+i] ■ ■ ■ Pk[ak\Xv 
= tr[Pfe-i+i[afe_;+i] • • • Pk[ak]Xvw'^] 



(26) 



for all . . . , Qf/g. Hence, if 

{Pk-i+i[(y^k-i+i] • '■Pk[(y^k]}cxk-i+u...,cxk spans 
(]-;;Dfc_i+ixDfc+i Q^g^ (j-;;^ ^.j^^g equivalent to Xvw^ = 0. 

Now, as It; 7^ (implied by tr[0] ^ 0), this is equivalent 
to Xv = 0. Hence, 

ker[r] = {X e I Xi; = 0} , (27) 

i.e., the rank of F is equal to 

^/c+i^/c+r+i - dim {X e (D^'^+i^^^'^+^+i \Xv = 0}. 

(28) 

Now, if {Pfe+l[Q^/c+l] • • •P/c+r[Q^fe+r]}afc + i,...,afc+, SpaUS 



ran[£;,^+;;-^+'-/]=ran[r], 
i.e., the rank of E^^'^^^'^'^'^^y is equal to 



(29) 



Dk+iDk+r+i - dim {X e (D^'=+ix-D'=+'-+i \Xv = 0}. 

(30) 
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As V ^ (implied by tr[0] ^ 0), we may set Vi = v 
and assume that there are vectors Vi G ([^^k+r+ixi ^ 
i = 2, . . . ,L)/e+r+i, such that {'yi}i=i,...,Dfc+,+i is an or- 
thogonal basis for (C^fc+r+ixi^ Letting {i/i}i=i,.. .,1)^+1 an 
orthogonal basis for (C^fc+ixi^ ^^y. 



X = ^ J2 XijUiv], 
i=i j=i 



(31) 



i.e., = Xv = Xvi is equivalent to = Xi^i for all 
z = 1, . . . , Dk-\-i. Hence, the rank of E^j^_i^'^' is equal 



to 



^fe+l^/c+r+l - ^fe+l(^fe+r+l — 1) = (32) 

Finally, 

rank[E{f+^-j-^>] <Dfc+i. (33) 



5. Now, assume that the errors of the local expectation 

values = (4"'^ • • • ^k+R-i^Q^ = 0' ^' ^' 

known for all k. Here, we simulate the error as a random 
variable X drawn from a Gaussian distribution with zero 
mean and standard deviation a. The errors in the mea- 
surements propagate into errors of the matrix A. In par- 
ticular, we find B = A + G, where Gx,y = X/Vd^ 
is a random matrix with entries drawn from a Gaus- 
sian distribution with zero mean and standard deviation 

Now, when applying the iterative scheme we have to 
solve linear equations of the form Bx = e where B is 
as described above and e is known. Instead of directly 
inverting this equation, we take possible variations in the 
matrix B into account, implying that the entries of B 
are themselves prone to noise and attempting to undo 
the imperfect measurements. This is done by solving the 
statistical least-squares problem [15] 



Non-invertible inputs 

The main issue arising when applying the iterative 
scheme to experimental data is that the local reductions 
are not known exactly. But of course, we may simply use 
their estimates (which are not necessarily states) as an in- 
put to compute the maps Ej^'^i'^~^^~^y and Ej^'^i'^~^^}_-^y. 
However, as we need to compute the inverse of the former 
map, already a small uncertainty will lead to a large er- 
ror in the inverse. This can be dealt with the method of 
stochastic robust approximation [15^, in particular with 
the Tikhonov regularization where a smoothing factor 
suppresses the effect of the smallest singular values of 
^{fc^.^.,/c4^r-^i} inverse. This is an often used method 

[16l[T7]. For completeness, we outline the derivation of 
the smoothing factor. Let us first ease notation a bit. 
We write 



k-l 



k-1 
)(afc+r-i) 



1 d?^ 
±, . . . , a , 



' X — 1, . . . , a . 



(34) 



Suppose now that one had access to the exact local 
expectation values. The matrix representation, A, of 

//j k-\-r 1| 

E}J"i' ,^ would then be given by 



^{k-l,...,k-l} 



. V f) 7-i{fe,...,/c+r— 1} / 

tT[PxQyd]. 



(35) 



Instead, we have only access to the noisy version of the 
entries tr[PxQyO]. Let us denote the resulting matrix by 



minE[||(B + G')a;-e||2] 



(36) 



where E denotes the expectation value and G' is a ma- 
trix with entries that are independent random numbers 
drawn from a Gaussian distribution with zero mean and 
standard deviation a j^f^^ ^ just like the assumed error 
in the measurements. This minimization problem can be 
rewritten as [T5] 



\Bx- 



\Vx\\ 



(37) 



with V = V'E[(G0^G']. Since we assume that 
the entries of G' are independent random numbers 
with zero mean and standard deviation a, we have 

problem reads 



<^a;,yC^ such that the minimization 



mm 



WBx- 



(38) 



Problems of this form are known as Tikhonov regulariza- 
tions which can be solved analytically [15, . The solution 
is given by 



X 



{B^B 



B'e. 



(39) 



Finally, let us denote the singular values of 5 by 5 1 > 
• • • > Srnin{d2^d2^} such that B = USV^ . Then, the solu- 



tion of the statistical least-square problem (36) is given 
by 5 = VSU^ where 5 is a diagonal matrix with entries 
fi/si where fi = s'^ / {s'^ + is the desired smoothing 
factor. 



